Animal dispersal costs are not universal

Supplmentary Material

Author

Martinig, A. R., S. L. P. Burk, Y. Yang, M. Lagisz, and S. Nakagawa

Published

February 10, 2026

1 Setting up

1.1 Load packages

Code
#force install pacman and orchaRd

#install.packages("pacman")
#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)
#devtools::install_github("daniel1noble/orchaRd", force = TRUE)

#load all other packages
pacman::p_load(tidyverse, # tidy family and related packages below
               kableExtra, # nice tables
               MuMIn,  # multi-model inference
               emmeans, # post-hoc tests
               gridExtra, # may not use this
               pander,   # nice tables
               metafor,  # package for meta-analysis
               ape,      # phylogenetic analysis
               ggtree,
               ggstance,
               MuMIn,  # multi-model inference
               patchwork,   # putting ggplots together - you need to install via devtool
               here,         # making reading files easy
               orchaRd, # plotting orchard plots
               matrixcalc, # matrix calculations
               cowplot #multipanel plots
)

eval(metafor:::.MuMIn)

1.2 Load data: paper and tree data

Code
#the dataset with effect sizes
dat <- read.csv(here("data", "clean_data.csv"))

#dim(dat)
#head(dat)

#the phylogenetic tree
tree <- read.tree(here("data", "species_tree.tre"))

1.3 Calculate effect sizes and sampling variances using custom functions

Code
# function to calculate effect sizes
# Zr - correlation
# there is always n

effect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments
                        est , se, p_val, direction, method){

  if(method == "mean_method"){
  
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r    
    
  }else if(method == "count_method"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- sqrt(m1)
    s2 <- sqrt(m2)
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r     
    
  }else if(method == "percent_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "percent_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    sd1 <- sd1/100
    sd2 <- sd2/100
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1))
    m2 <- asin(sqrt(m2))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "t_method1"){
    

    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "t_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "F_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
  
  }else if(method == "F_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "correlation_method1"){
    
    r <- est
    
  }else if(method == "correlation_method2"){
    
    r <- 2*sin((pi/6)*est)
    
  }else if(method == "estimate_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  }else if(method == "estimate_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  } 
  
    if(r >= 1){
    # if over 1, we use 0.99
    Zr <- atanh(0.99)

    }else if(r <= -1){

    Zr <- atanh(-0.99) # r = correlation

    } else {

    Zr <- atanh(r) # r = correlation

    }
  
  VZr <- 1 /(n - 3)
  
  data.frame(ri = r, yi = Zr , vi = VZr)
  
}


##########

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name","\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "normal") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  # creating a function to run models
  run_rma1 <- function(name) {
      
    # variance covarinace matrix for sampling errors
    VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
    #VCV1 <- nearPD(VCV1)$mat
      
    rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }

    run_rma2 <- function(name) {
    
            # variance covarinace matrix for sampling errors
            VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
            #VCVa <- nearPD(VCVa)$mat
               
               rma.mv(abs_yi2, V = VCVa,
               mods = ~relevel(dat2[[mod]], ref = name),
                      random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
    }
    
    run_rma3 <- function(name) {
      
              # variance covarinace matrix for sampling errors
                VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
                #VCV1 <- nearPD(VCV1)$mat
               
                rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list( # hetero in relation to mod
                            formula(paste("~", mod, "| effectID")),
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     struct = "DIAG",
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }


# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
# use different specifications for aboslute and hetero models 
if (type == "normal"){

    model_all <- purrr::map(level_names[-length(level_names)], run_rma1)

  }else if(type == "absolute"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma2)
    
  }else if(type == "hetero"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma3)
    
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# absolute value analyses

# folded mean
folded_mu <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_mu
} 

# folded variance
folded_v <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
  # adding se to make bigger mean
  fold_v <-fold_se^2
  fold_v
} 

1.4 Preparing final dataset

2 Intercept meta-analytic model

2.1 Phylogenetic multilevel models

Code
# takes a while to run
mod <- rma.mv(yi = yi, 
              V = VCV,
              mod = ~ 1,
              data = dat,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned,
                            ~ 1 | phylogeny),
              R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)

# saving the runs
saveRDS(mod, here("Rdata", "mod.RDS"))
Code
# getting the saved model
mod <- readRDS(here("Rdata", "mod.RDS"))

summary(mod)

Multivariate Meta-Analysis Model (k = 715; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.1889   486.3779   496.3779   519.2323   496.4626   

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R 
sigma^2.1  0.0734  0.2709    715     no         effectID   no 
sigma^2.2  0.0207  0.1438    209     no          paperID   no 
sigma^2.3  0.0175  0.1323    149     no  species_cleaned   no 
sigma^2.4  0.0000  0.0000    149     no        phylogeny  yes 

Test for Heterogeneity:
Q(df = 714) = 15803.8813, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0244  0.0294  -0.8311  714  0.4062  -0.0822  0.0333    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# Phylogeney accounts for nothing (0.00%), so we can take it out. 
round(i2_ml(mod),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.01              63.81              17.98              15.22 
      I2_phylogeny 
              0.00 
Code
pred_mod <- predict.rma(mod)

  estimate <- pred_mod$pred
  lowerCL <- pred_mod$ci.lb
  upperCL <- pred_mod$ci.ub 
  lowerPR <- pred_mod$cr.lb
  upperPR <- pred_mod$cr.ub 
  
  table_mod <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod <- table_mod %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod as RDS
saveRDS(tab_mod, here("Rdata", "tab_mod.RDS"))

Table S1. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for intercept meta-analytic model with all five random effects.

Code
tab_mod <-readRDS(here("Rdata", "tab_mod.RDS"))

tab_mod
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.024 -0.082 0.033 0.406 -0.683 0.634

We remove phylogeny as a random effect from our meta-analytic model because phylogeny accounts for nothing (0%)

Code
mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.3245   486.6490   494.6490   512.9437   494.7053   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0733  0.2708    717     no         effectID 
sigma^2.2  0.0210  0.1448    210     no          paperID 
sigma^2.3  0.0167  0.1294    150     no  species_cleaned 

Test for Heterogeneity:
Q(df = 716) = 15803.9745, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0241  0.0293  -0.8249  716  0.4097  -0.0816  0.0333    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod1),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             96.99              64.05              18.32              14.62 
Code
reduced_intercept<-orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)
reduced_intercept

Figure S1. Overall meta-analytic mean effect size (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE).
Code
pred_mod1 <- predict.rma(mod1)

  estimate <- pred_mod1$pred
  lowerCL <- pred_mod1$ci.lb
  upperCL <- pred_mod1$ci.ub 
  lowerPR <- pred_mod1$cr.lb
  upperPR <- pred_mod1$cr.ub 
  
  table_mod1 <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod1$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod1 <- table_mod1 %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod1 as RDS
saveRDS(tab_mod1, here("Rdata", "tab_mod1.RDS"))

Table S2. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for reduced intercept meta-analytic model without phylogeny as a random effect.

Code
tab_mod1 <-readRDS(here("Rdata", "tab_mod1.RDS"))

tab_mod1
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.024 -0.082 0.033 0.41 -0.681 0.633

3 Uni-moderator meta-regression models

For each of our moderators, we run a uni-moderator meta-regression model

Code
mod_class <- rma.mv(yi = yi, V = VCV,
               mod = ~ species_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_class)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.0715   482.1429   500.1429   541.2430   500.3997   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0724  0.2691    717     no         effectID 
sigma^2.2  0.0151  0.1227    210     no          paperID 
sigma^2.3  0.0293  0.1712    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 711) = 15272.7596, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 711) = 0.5852, p-val = 0.7423

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0222  0.1645   0.1348  711  0.8928  -0.3009 
species_classArachnida         0.0412  0.2679   0.1537  711  0.8779  -0.4847 
species_classAves             -0.0288  0.0356  -0.8085  711  0.4191  -0.0986 
species_classInsecta           0.0955  0.1379   0.6920  711  0.4891  -0.1754 
species_classMammalia         -0.0560  0.0444  -1.2600  711  0.2081  -0.1431 
species_classReptilia          0.0938  0.1204   0.7788  711  0.4363  -0.1426 
                              ci.ub    
species_classActinopterygii  0.3452    
species_classArachnida       0.5671    
species_classAves            0.0411    
species_classInsecta         0.3663    
species_classMammalia        0.0312    
species_classReptilia        0.3302    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_class)
   R2_marginal R2_conditional 
   0.008721663    0.385329063 
Code
orchard_plot(mod_class, mod = "species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S2. Mean fitness effect of different taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_class <- all_models(mod_class,  mod = "species_class", type = "normal")

# saving tab_class as RDS
saveRDS(tab_class, here("Rdata", "tab_class.RDS"))

Table S3. Mean estimates of different taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_class <-readRDS(here("Rdata", "tab_class.RDS"))

tab_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Actinopterygii 0.022 -0.302 0.347 0.893 -0.725 0.769
Arachnida 0.041 -0.487 0.569 0.879 -0.814 0.896
Aves -0.029 -0.099 0.041 0.419 -0.705 0.647
Insecta 0.096 -0.175 0.368 0.487 -0.629 0.822
Mammalia -0.056 -0.144 0.031 0.207 -0.735 0.622
Reptilia 0.094 -0.144 0.331 0.440 -0.620 0.807
Actinopterygii-Arachnida 0.019 -0.598 0.636 0.952 -0.894 0.932
Actinopterygii-Aves -0.051 -0.378 0.276 0.759 -0.799 0.697
Actinopterygii-Insecta 0.074 -0.345 0.493 0.729 -0.718 0.866
Actinopterygii-Mammalia -0.079 -0.408 0.250 0.639 -0.827 0.670
Actinopterygii-Reptilia 0.071 -0.324 0.467 0.723 -0.709 0.852
Arachnida-Aves -0.070 -0.600 0.460 0.796 -0.927 0.787
Arachnida-Insecta 0.055 -0.536 0.646 0.855 -0.841 0.951
Arachnida-Mammalia -0.098 -0.629 0.434 0.719 -0.955 0.760
Arachnida-Reptilia 0.053 -0.523 0.628 0.858 -0.833 0.938
Aves-Insecta 0.125 -0.151 0.401 0.373 -0.602 0.852
Aves-Mammalia -0.028 -0.123 0.068 0.570 -0.707 0.652
Aves-Reptilia 0.123 -0.118 0.363 0.317 -0.592 0.837
Insecta-Mammalia -0.153 -0.431 0.126 0.282 -0.881 0.575
Insecta-Reptilia -0.002 -0.358 0.353 0.989 -0.763 0.758
Mammalia-Reptilia 0.150 -0.093 0.393 0.225 -0.565 0.865
Code
mod_type <- rma.mv(yi = yi,  
               V = VCV,
               mod = ~ study_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_type)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0499   486.0999   496.0999   518.9613   496.1845   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0732  0.2705    717     no         effectID 
sigma^2.2  0.0202  0.1421    210     no          paperID 
sigma^2.3  0.0187  0.1366    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15759.6360, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 715) = 0.3938, p-val = 0.6746

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0236  0.0300  -0.7879  715  0.4310  -0.0825  0.0352 
study_typesemi-natural   -0.0408  0.0659  -0.6195  715  0.5358  -0.1703  0.0886 
                          
study_typenatural         
study_typesemi-natural    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_type)
   R2_marginal R2_conditional 
  0.0002335152   0.3469633182 
Code
orchard_plot(mod_type, mod = "study_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S3. Mean fitness effect of natural and semi-natural study types. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_type <- all_models(mod_type,  mod = "study_type", type = "normal")

# saving tab_type as RDS
saveRDS(tab_type, here("Rdata", "tab_type.RDS"))

Table S4. Mean estimates of natural and semi-natural study types with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_type <-readRDS(here("Rdata", "tab_type.RDS"))

tab_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Natural -0.024 -0.083 0.035 0.428 -0.685 0.638
Semi-natural -0.042 -0.171 0.088 0.529 -0.713 0.630
Natural-Semi-natural -0.018 -0.144 0.108 0.781 -0.689 0.653
Code
mod_design <- rma.mv(yi = yi, V = VCV,
               mod = ~ study_design - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_design)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0290   486.0580   496.0580   518.9194   496.1426   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0733  0.2707    717     no         effectID 
sigma^2.2  0.0232  0.1523    210     no          paperID 
sigma^2.3  0.0144  0.1201    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15788.5009, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 715) = 0.4376, p-val = 0.6457

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast   -0.0192  0.0303  -0.6339  715  0.5264  -0.0788  0.0403 
study_designgradient   -0.0398  0.0448  -0.8877  715  0.3750  -0.1277  0.0482 
                        
study_designcontrast    
study_designgradient    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_design)
   R2_marginal R2_conditional 
  0.0006775013   0.3396059825 
Code
orchard_plot(mod_design, mod = "study_design", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S4. Mean fitness effect of gradient and contrast study designs. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_design <- all_models(mod_design,  mod = "study_design", type = "normal")

# saving tab_design as RDS
saveRDS(tab_design, here("Rdata", "tab_design.RDS"))

Table S5. Mean estimates of gradient and contrast study designs with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_design <-readRDS(here("Rdata", "tab_design.RDS"))

tab_design
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Contrast -0.02 -0.079 0.040 0.520 -0.678 0.638
Gradient -0.04 -0.128 0.048 0.374 -0.701 0.621
Contrast-Gradient -0.02 -0.104 0.063 0.633 -0.681 0.640
Code
mod_disp <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ dispersal_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.8184   485.6368   497.6368   525.0621   497.7556   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0733  0.2708    717     no         effectID 
sigma^2.2  0.0225  0.1501    210     no          paperID 
sigma^2.3  0.0155  0.1244    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 714) = 15631.5675, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 714) = 0.4324, p-val = 0.7299

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeBoth       -0.0812  0.0795  -1.0219  714  0.3072  -0.2372  0.0748 
dispersal_typebreeding   -0.0142  0.0418  -0.3389  714  0.7348  -0.0962  0.0679 
dispersal_typenatal      -0.0221  0.0317  -0.6975  714  0.4857  -0.0843  0.0401 
                          
dispersal_typeBoth        
dispersal_typebreeding    
dispersal_typenatal       

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp)
   R2_marginal R2_conditional 
   0.001371961    0.342378689 
Code
orchard_plot(mod_disp, mod = "dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S5. Mean fitness effect of natal and breeding dispersal (or both dispersal types grouped). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_disp <- all_models(mod_disp,  mod = "dispersal_type", type = "normal")

# saving tab_disp as RDS
saveRDS(tab_disp, here("Rdata", "tab_disp.RDS"))

Table S6. Mean estimates of natal and breeding dispersal (or both dispersal types combined) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_disp <-readRDS(here("Rdata", "tab_disp.RDS"))

tab_disp
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Both -0.081 -0.238 0.075 0.307 -0.756 0.594
Breeding -0.014 -0.097 0.068 0.731 -0.676 0.647
Natal -0.022 -0.085 0.040 0.482 -0.682 0.637
Both-Breeding 0.067 -0.100 0.234 0.431 -0.610 0.744
Both-Natal 0.059 -0.097 0.215 0.458 -0.616 0.734
Breeding-Natal -0.008 -0.088 0.072 0.845 -0.669 0.653
Code
mod_phase <- rma.mv(yi = yi, V = VCV,
               mod = ~ dispersal_phase - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_phase)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.5522   485.1043   495.1043   517.9657   495.1889   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0731  0.2703    717     no         effectID 
sigma^2.2  0.0217  0.1474    210     no          paperID 
sigma^2.3  0.0164  0.1279    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15801.2403, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 715) = 1.0664, p-val = 0.3448

Model Results:

                            estimate      se     tval   df    pval    ci.lb 
dispersal_phaseprospection   -0.0632  0.0437  -1.4470  715  0.1483  -0.1490 
dispersal_phasesettlement    -0.0145  0.0304  -0.4770  715  0.6335  -0.0741 
                             ci.ub    
dispersal_phaseprospection  0.0226    
dispersal_phasesettlement   0.0451    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_phase)
   R2_marginal R2_conditional 
   0.003188043    0.344646396 
Code
orchard_plot(mod_phase, mod = "dispersal_phase", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S6. Mean fitness effect of prospection and settlement dispersal phases. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_phase <- all_models(mod_phase,  mod = "dispersal_phase", type = "normal")

# saving tab_phase as RDS
saveRDS(tab_phase, here("Rdata", "tab_phase.RDS"))

Table S7. Mean estimates of prospection and settlement dispersal phases with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_phase <-readRDS(here("Rdata", "tab_phase.RDS"))

tab_phase
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Prospection -0.067 -0.158 0.025 0.152 -0.836 0.703
Settlement -0.022 -0.085 0.041 0.493 -0.788 0.745
Prospection-Settlement 0.045 -0.042 0.132 0.314 -0.724 0.813
Code
mod_sex <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.8528   483.7055   495.7055   523.1308   495.8243   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0726  0.2694    717     no         effectID 
sigma^2.2  0.0260  0.1613    210     no          paperID 
sigma^2.3  0.0118  0.1087    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 714) = 15694.0840, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 714) = 1.3896, p-val = 0.2447

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
sexBoth     -0.0770  0.0410  -1.8779  714  0.0608  -0.1574  0.0035  . 
sexFemale   -0.0056  0.0327  -0.1708  714  0.8644  -0.0698  0.0587    
sexMale     -0.0002  0.0337  -0.0071  714  0.9944  -0.0665  0.0660    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex)
   R2_marginal R2_conditional 
   0.007006424    0.347383516 
Code
orchard_plot(mod_sex, mod = "sex", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S7. Mean fitness effect of females, males, or both sexes grouped. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex <- all_models(mod_sex,  mod = "sex", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))

Table S8. Mean estimates of females, males, or both sexes grouped with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))

tab_sex
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.073 -0.162 0.015 0.102 -0.841 0.694
F -0.013 -0.083 0.056 0.711 -0.779 0.753
M -0.017 -0.088 0.055 0.650 -0.782 0.749
B-F 0.060 -0.032 0.152 0.199 -0.708 0.828
B-M 0.057 -0.037 0.151 0.237 -0.711 0.825
F-M -0.003 -0.062 0.055 0.908 -0.768 0.761

3.1 Broad level

Code
mod_age1 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age1)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.7487   483.4974   495.4974   522.9227   495.6162   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0736  0.2714    717     no         effectID 
sigma^2.2  0.0173  0.1317    210     no          paperID 
sigma^2.3  0.0204  0.1430    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 714) = 15749.9440, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 714) = 1.3708, p-val = 0.2505

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0196  0.0304  -0.6464  714  0.5182  -0.0793 
age_class_cleanjuvenile    0.0065  0.0507   0.1283  714  0.8979  -0.0929 
age_class_cleanmix        -0.1139  0.0588  -1.9355  714  0.0533  -0.2294 
                          ci.ub    
age_class_cleanadult     0.0400    
age_class_cleanjuvenile  0.1059    
age_class_cleanmix       0.0016  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age1)
   R2_marginal R2_conditional 
   0.005819041    0.342929616 
Code
orchard_plot(mod_age1, mod = "age_class_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S8. Mean fitness effect of life history stage (adult, juvenile, or mixed). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_age1 <- all_models(mod_age1,  mod = "age_class_clean", type = "normal")

# saving tab_age1 as RDS
saveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))

Table S9. Mean estimates of life history stage (adult, juvenile, or mixed) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))

tab_age1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Adult -0.004 -0.061 0.052 0.879 -0.735 0.726
Juvenile 0.012 -0.080 0.105 0.794 -0.722 0.747
Mix -0.101 -0.224 0.022 0.108 -0.839 0.638
Adult-Juvenile 0.017 -0.072 0.106 0.713 -0.717 0.751
Adult-Mix -0.096 -0.217 0.025 0.118 -0.835 0.642
Juvenile-Mix -0.113 -0.254 0.027 0.115 -0.855 0.629

3.2 Finer level

Code
# age_class

mod_age2 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age2)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.5957   481.1915   501.1915   546.8441   501.5062   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0741  0.2721    717     no         effectID 
sigma^2.2  0.0188  0.1371    210     no          paperID 
sigma^2.3  0.0180  0.1340    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 710) = 15553.3665, p-val < .0001

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 710) = 0.8752, p-val = 0.5257

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
age_classA     -0.0115  0.0330  -0.3486  710  0.7275  -0.0763  0.0533    
age_classJ      0.0053  0.0508   0.1040  710  0.9172  -0.0945  0.1050    
age_classJA    -0.3023  0.1781  -1.6979  710  0.0900  -0.6519  0.0473  . 
age_classJY    -0.1218  0.0834  -1.4605  710  0.1446  -0.2856  0.0419    
age_classJYA   -0.0654  0.0846  -0.7726  710  0.4400  -0.2315  0.1007    
age_classY     -0.0439  0.0567  -0.7734  710  0.4396  -0.1552  0.0675    
age_classYA    -0.0374  0.0460  -0.8113  710  0.4175  -0.1278  0.0530    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age2)
   R2_marginal R2_conditional 
     0.0110822      0.3391335 
Code
orchard_plot(mod_age2, mod = "age_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S9. Mean fitness effect of life history stage at the level of precision reported (A: adult, Y: yearling, J: juvenile). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)

3.3 Broad level

Code
mod_fit1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_higher_level - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit1)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.0456   484.0912   496.0912   523.5165   496.2100   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0731  0.2703    717     no         effectID 
sigma^2.2  0.0209  0.1445    210     no          paperID 
sigma^2.3  0.0169  0.1300    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 714) = 15677.2739, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 714) = 1.3311, p-val = 0.2631

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2390  0.1437  -1.6630  714  0.0968 
fitness_higher_levelreproduction       -0.0092  0.0317  -0.2893  714  0.7724 
fitness_higher_levelsurvival           -0.0353  0.0324  -1.0892  714  0.2764 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5212  0.0432  . 
fitness_higher_levelreproduction      -0.0715  0.0531    
fitness_higher_levelsurvival          -0.0989  0.0283    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit1)
   R2_marginal R2_conditional 
    0.00446575     0.34368738 
Code
orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S10. Mean fitness effect of different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_fit1 <- all_models(mod_fit1,  mod = "fitness_higher_level", type = "normal")

# saving tab_fit1 as RDS
saveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))

Table S10. Mean estimates of different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))

tab_fit1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Indirect fitness -0.237 -0.572 0.098 0.166 -1.068 0.595
Reproduction -0.019 -0.084 0.047 0.580 -0.782 0.745
Survival -0.038 -0.106 0.030 0.273 -0.802 0.726
Indirect fitness-Reproduction 0.218 -0.117 0.554 0.202 -0.613 1.050
Indirect fitness-Survival 0.199 -0.136 0.534 0.244 -0.632 1.030
Reproduction-Survival -0.019 -0.079 0.040 0.525 -0.783 0.744

3.4 Finer level

Code
mod_fit2 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_metric_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit2)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.8051   481.6102   503.6102   553.8126   503.9890   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0723  0.2689    717     no         effectID 
sigma^2.2  0.0163  0.1277    210     no          paperID 
sigma^2.3  0.0254  0.1593    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 709) = 15212.0108, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 709) = 1.0802, p-val = 0.3749

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0822  0.0855 
fitness_metric_cleanindirect fitness                        -0.2326  0.1443 
fitness_metric_cleanlifetime breeding success                0.0070  0.0367 
fitness_metric_cleanlifetime reproductive success           -0.0258  0.0393 
fitness_metric_cleanoffspring reproduction                   0.0345  0.0867 
fitness_metric_cleanoffspring survival                       0.0076  0.0437 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0207  0.0836 
fitness_metric_cleansurvival                                -0.0602  0.0361 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.9610  709  0.3369 
fitness_metric_cleanindirect fitness                       -1.6118  709  0.1074 
fitness_metric_cleanlifetime breeding success               0.1905  709  0.8490 
fitness_metric_cleanlifetime reproductive success          -0.6567  709  0.5116 
fitness_metric_cleanoffspring reproduction                  0.3973  709  0.6913 
fitness_metric_cleanoffspring survival                      0.1749  709  0.8612 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.2475  709  0.8046 
fitness_metric_cleansurvival                               -1.6667  709  0.0960 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2501  0.0857    
fitness_metric_cleanindirect fitness                       -0.5160  0.0507    
fitness_metric_cleanlifetime breeding success              -0.0650  0.0790    
fitness_metric_cleanlifetime reproductive success          -0.1030  0.0514    
fitness_metric_cleanoffspring reproduction                 -0.1358  0.2048    
fitness_metric_cleanoffspring survival                     -0.0781  0.0934    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.1848  0.1434    
fitness_metric_cleansurvival                               -0.1311  0.0107  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit2)
   R2_marginal R2_conditional 
    0.01110422     0.37273928 
Code
orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 0)

Figure S11. Mean fitness effect of different fitness components at a finer resolution. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)
Code
mod_gen <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ whose_fitness - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_gen)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.8208   485.6416   495.6416   518.5031   495.7263   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0728  0.2698    717     no         effectID 
sigma^2.2  0.0179  0.1337    210     no          paperID 
sigma^2.3  0.0224  0.1495    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15803.7544, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 715) = 1.0797, p-val = 0.3403

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessdescendant    0.0085  0.0421   0.2024  715  0.8397  -0.0741 
whose_fitnessfocal        -0.0320  0.0303  -1.0567  715  0.2910  -0.0915 
                          ci.ub    
whose_fitnessdescendant  0.0912    
whose_fitnessfocal       0.0275    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_gen)
   R2_marginal R2_conditional 
    0.00198268     0.35732506 
Code
orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S12. Mean fitness effect of focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_gen <- all_models(mod_gen,  mod = "whose_fitness", type = "normal")

# saving tab_gen as RDS
saveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))

Table S11. Mean estimates of non-dispersers and dispersers and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))

tab_gen
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Individual -0.038 -0.100 0.023 0.219 -0.808 0.731
Offspring 0.005 -0.086 0.096 0.907 -0.767 0.777
Individual-Offspring 0.044 -0.035 0.123 0.275 -0.727 0.815

3.5 Normal analysis

Code
# getting absolute values for effect sizes as we expect shorter years will have more fluctuations
dat$ln_study_duration <- log(dat$study_duration+1)

mod_dur <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ ln_study_duration,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0516   486.1032   496.1032   518.9646   496.1878   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0733  0.2707    717     no         effectID 
sigma^2.2  0.0210  0.1448    210     no          paperID 
sigma^2.3  0.0174  0.1320    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15752.6262, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 715) = 0.0185, p-val = 0.8917

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0176  0.0584  -0.3018  715  0.7629  -0.1323  0.0970    
ln_study_duration   -0.0028  0.0205  -0.1362  715  0.8917  -0.0431  0.0375    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_dur)
   R2_marginal R2_conditional 
  6.279703e-05   3.438337e-01 
Code
bubble_plot(mod_dur,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S13. Mean fitness effect of study duration (log years). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# absolute value analysis
mod_dur_abs <- rma.mv(yi = yi_abs,
                V = VCVa,
               mod = ~ ln_study_duration,
               data = dat,
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur_abs)

Multivariate Meta-Analysis Model (k = 717; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
 38.0638  -76.1275  -66.1275  -43.2661  -66.0429   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0316  0.1779    717     no         effectID 
sigma^2.2  0.0070  0.0837    210     no          paperID 
sigma^2.3  0.0107  0.1036    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 9300.5284, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 715) = 3.5428, p-val = 0.0602

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0336  0.0390  -0.8600  715  0.3901  -0.1102  0.0431    
ln_study_duration    0.0256  0.0136   1.8822  715  0.0602  -0.0011  0.0523  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_dur_abs)
   R2_marginal R2_conditional 
    0.01178184     0.36672938 
Code
bubble_plot(mod_dur_abs,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S14. Mean fitness effect of study duration (log years) from absolute value analysis. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

3.6 Normal analysis

Code
mod_focus <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focus)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0401   486.0802   496.0802   518.9416   496.1648   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0731  0.2704    717     no         effectID 
sigma^2.2  0.0209  0.1446    210     no          paperID 
sigma^2.3  0.0178  0.1335    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15784.2497, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 715) = 0.4505, p-val = 0.6375

Model Results:

                                 estimate      se     tval   df    pval 
fitness_main_focusNon-selective   -0.0095  0.0451  -0.2103  715  0.8335 
fitness_main_focusSelective       -0.0285  0.0307  -0.9289  715  0.3533 
                                   ci.lb   ci.ub    
fitness_main_focusNon-selective  -0.0981  0.0791    
fitness_main_focusSelective      -0.0887  0.0317    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus)
   R2_marginal R2_conditional 
  0.0005755634   0.3466375060 
Code
orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S15. Mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_focus <- all_models(mod_focus,  mod = "fitness_main_focus", type = "normal")

# saving tab_focus as RDS
saveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))

Table S12. Mean estimates of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))

tab_focus
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
N -0.019 -0.110 0.073 0.690 -0.787 0.750
Y -0.034 -0.096 0.029 0.292 -0.799 0.732
N-Y -0.015 -0.103 0.073 0.734 -0.783 0.752

3.7 Comparing normal model to heteroscedastic model

Code
# homo
mod_focus1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_focus1)

Multivariate Meta-Analysis Model (k = 717; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-243.9040  2928.1024   497.8080   520.6834   497.8924   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0733  0.2707    717     no         effectID 
sigma^2.2  0.0211  0.1454    210     no          paperID 
sigma^2.3  0.0162  0.1273    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15784.2497, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 715) = 0.4291, p-val = 0.6512

Model Results:

                                 estimate      se     tval   df    pval 
fitness_main_focusNon-selective   -0.0092  0.0448  -0.2052  715  0.8375 
fitness_main_focusSelective       -0.0276  0.0304  -0.9071  715  0.3646 
                                   ci.lb   ci.ub    
fitness_main_focusNon-selective  -0.0972  0.0788    
fitness_main_focusSelective      -0.0872  0.0321    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus1)
   R2_marginal R2_conditional 
  0.0005423827   0.3379895506 
Code
# hetero
mod_focus2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ fitness_main_focus,
                 data = dat, 
                 random = list(
                   ~ fitness_main_focus | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_focus2)

Multivariate Meta-Analysis Model (k = 717; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-239.0021  2918.2986   490.0042   517.4546   490.1225   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0234  0.1528    210     no          paperID 
sigma^2.2  0.0158  0.1258    150     no  species_cleaned 

outer factor: effectID           (nlvls = 717)
inner factor: fitness_main_focus (nlvls = 2)

            estim    sqrt  k.lvl  fixed          level 
tau^2.1    0.0425  0.2063    166     no  Non-selective 
tau^2.2    0.0800  0.2828    551     no      Selective 

Test for Residual Heterogeneity:
QE(df = 715) = 15784.2497, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 715) = 0.2523, p-val = 0.6156

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                       -0.0054  0.0423  -0.1271  715  0.8989  -0.0885 
fitness_main_focusSelective   -0.0203  0.0405  -0.5023  715  0.6156  -0.0999 
                              ci.ub    
intrcpt                      0.0778    
fitness_main_focusSelective  0.0592    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
#r2_ml(mod_focusA1)

# they are not significantly different
anova(mod_focus1, mod_focus2)

        df      AIC      BIC     AICc    logLik    LRT   pval         QE 
Full     6 490.0042 517.4546 490.1225 -239.0021               15784.2497 
Reduced  5 497.8080 520.6834 497.8924 -243.9040 9.8039 0.0017 15784.2497 
Code
# homo
orchard_plot(mod_focus1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S16. Maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# hetero
orchard_plot(mod_focus2, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S17. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

3.8 Normal analysis

Code
mod_confirm <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_confirm)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.3569   484.7138   498.7138   530.7002   498.8727   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0731  0.2703    717     no         effectID 
sigma^2.2  0.0275  0.1658    210     no          paperID 
sigma^2.3  0.0107  0.1032    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 713) = 15676.2331, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 713) = 0.3593, p-val = 0.8376

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0294  0.0819   0.3590  713  0.7197  -0.1315 
confirmation_biasCost      -0.0373  0.0383  -0.9724  713  0.3312  -0.1125 
confirmation_biasMixed     -0.0227  0.0381  -0.5960  713  0.5514  -0.0976 
confirmation_biasNeutral   -0.0038  0.0458  -0.0841  713  0.9330  -0.0938 
                           ci.ub    
confirmation_biasBenefit  0.1903    
confirmation_biasCost     0.0380    
confirmation_biasMixed    0.0521    
confirmation_biasNeutral  0.0861    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm)
   R2_marginal R2_conditional 
   0.002373472    0.344489760 
Code
orchard_plot(mod_confirm, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S18. Mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_confirm <- all_models(mod_confirm,  mod = "confirmation_bias", type = "normal")

# saving tab_confirm as RDS
saveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))

Table S13. Mean estimates of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_confirm <-readRDS(here("Rdata", "tab_confirm.RDS"))

tab_confirm
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Benefit 0.012 -0.154 0.178 0.884 -0.764 0.788
Cost -0.042 -0.122 0.037 0.296 -0.804 0.720
Mixed -0.028 -0.107 0.050 0.483 -0.790 0.734
Neutral -0.009 -0.102 0.084 0.852 -0.772 0.755
Benefit-Cost -0.055 -0.227 0.118 0.536 -0.832 0.723
Benefit-Mixed -0.040 -0.212 0.132 0.645 -0.818 0.737
Benefit-Neutral -0.021 -0.198 0.156 0.815 -0.800 0.757
Cost-Mixed 0.014 -0.081 0.110 0.770 -0.750 0.778
Cost-Neutral 0.033 -0.071 0.138 0.531 -0.732 0.799
Mixed-Neutral 0.019 -0.085 0.124 0.719 -0.746 0.784

3.9 Comparing normal model to heteroscedastic model

Code
# homo
mod_confirm1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_confirm1)

Multivariate Meta-Analysis Model (k = 717; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-243.5898  2927.4740   501.1796   533.2051   501.3376   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0733  0.2708    717     no         effectID 
sigma^2.2  0.0266  0.1632    210     no          paperID 
sigma^2.3  0.0092  0.0961    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 713) = 15676.2331, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 713) = 0.3488, p-val = 0.8449

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0294  0.0809   0.3638  713  0.7161  -0.1295 
confirmation_biasCost      -0.0363  0.0378  -0.9598  713  0.3375  -0.1104 
confirmation_biasMixed     -0.0211  0.0374  -0.5637  713  0.5731  -0.0946 
confirmation_biasNeutral   -0.0031  0.0452  -0.0678  713  0.9460  -0.0917 
                           ci.ub    
confirmation_biasBenefit  0.1884    
confirmation_biasCost     0.0379    
confirmation_biasMixed    0.0524    
confirmation_biasNeutral  0.0856    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm1)
   R2_marginal R2_conditional 
   0.002350484    0.329906807 
Code
# hetero
mod_confirm2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ confirmation_bias,
                 data = dat, 
                 random = list(
                   ~ confirmation_bias | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_confirm2)

Multivariate Meta-Analysis Model (k = 717; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-224.9869  2890.2683   469.9738   515.7246   470.2855   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0283  0.1681    210     no          paperID 
sigma^2.2  0.0031  0.0559    150     no  species_cleaned 

outer factor: effectID          (nlvls = 717)
inner factor: confirmation_bias (nlvls = 4)

            estim    sqrt  k.lvl  fixed    level 
tau^2.1    0.0130  0.1140     32     no  Benefit 
tau^2.2    0.0791  0.2813    210     no     Cost 
tau^2.3    0.0553  0.2351    297     no    Mixed 
tau^2.4    0.1295  0.3598    178     no  Neutral 

Test for Residual Heterogeneity:
QE(df = 713) = 15676.2331, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 713) = 0.6269, p-val = 0.5978

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
intrcpt                     0.0509  0.0625   0.8149  713  0.4154  -0.0718 
confirmation_biasCost      -0.0904  0.0665  -1.3588  713  0.1746  -0.2209 
confirmation_biasMixed     -0.0734  0.0653  -1.1241  713  0.2614  -0.2017 
confirmation_biasNeutral   -0.0634  0.0714  -0.8883  713  0.3747  -0.2035 
                           ci.ub    
intrcpt                   0.1737    
confirmation_biasCost     0.0402    
confirmation_biasMixed    0.0548    
confirmation_biasNeutral  0.0767    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# they are not significantly different
anova(mod_confirm1, mod_confirm2)

        df      AIC      BIC     AICc    logLik     LRT   pval         QE 
Full    10 469.9738 515.7246 470.2855 -224.9869                15676.2331 
Reduced  7 501.1796 533.2051 501.3376 -243.5898 37.2058 <.0001 15676.2331 
Code
# homo
orchard_plot(mod_confirm1, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S19. Maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# hetero
orchard_plot(mod_confirm2, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S20. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

4 Interaction meta-regression models (2 moderators)

Code
# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach 

dat$sex_species_class <- as.factor(paste(dat$sex, dat$species_class ,sep = "_"))

mod_sex_species_class <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_species_class - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_species_class)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-236.9614   473.9228   509.9228   591.8936   510.9242   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0727  0.2697    717     no         effectID 
sigma^2.2  0.0199  0.1412    210     no          paperID 
sigma^2.3  0.0194  0.1392    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 702) = 14581.0769, p-val < .0001

Test of Moderators (coefficients 1:15):
F(df1 = 15, df2 = 702) = 0.9697, p-val = 0.4858

Model Results:

                                        estimate      se     tval   df    pval 
sex_species_classBoth_Aves               -0.0231  0.0505  -0.4579  702  0.6472 
sex_species_classBoth_Insecta            -0.1449  0.2378  -0.6093  702  0.5425 
sex_species_classBoth_Mammalia           -0.1782  0.0639  -2.7880  702  0.0054 
sex_species_classBoth_Reptilia            0.1524  0.2705   0.5632  702  0.5735 
sex_species_classFemale_Actinopterygii   -0.0609  0.1926  -0.3163  702  0.7519 
sex_species_classFemale_Arachnida         0.0398  0.2585   0.1539  702  0.8777 
sex_species_classFemale_Aves             -0.0293  0.0405  -0.7246  702  0.4690 
sex_species_classFemale_Insecta           0.2478  0.1746   1.4190  702  0.1563 
sex_species_classFemale_Mammalia         -0.0146  0.0489  -0.2989  702  0.7651 
sex_species_classFemale_Reptilia          0.0864  0.1277   0.6766  702  0.4989 
sex_species_classMale_Actinopterygii      0.0937  0.1882   0.4980  702  0.6186 
sex_species_classMale_Aves               -0.0330  0.0407  -0.8113  702  0.4174 
sex_species_classMale_Insecta            -0.1122  0.3322  -0.3377  702  0.7357 
sex_species_classMale_Mammalia            0.0201  0.0534   0.3771  702  0.7062 
sex_species_classMale_Reptilia            0.0772  0.1577   0.4895  702  0.6246 
                                          ci.lb    ci.ub     
sex_species_classBoth_Aves              -0.1224   0.0761     
sex_species_classBoth_Insecta           -0.6117   0.3219     
sex_species_classBoth_Mammalia          -0.3036  -0.0527  ** 
sex_species_classBoth_Reptilia          -0.3788   0.6835     
sex_species_classFemale_Actinopterygii  -0.4389   0.3171     
sex_species_classFemale_Arachnida       -0.4677   0.5473     
sex_species_classFemale_Aves            -0.1088   0.0502     
sex_species_classFemale_Insecta         -0.0951   0.5906     
sex_species_classFemale_Mammalia        -0.1106   0.0814     
sex_species_classFemale_Reptilia        -0.1643   0.3371     
sex_species_classMale_Actinopterygii    -0.2758   0.4633     
sex_species_classMale_Aves              -0.1130   0.0469     
sex_species_classMale_Insecta           -0.7644   0.5400     
sex_species_classMale_Mammalia          -0.0847   0.1250     
sex_species_classMale_Reptilia          -0.2324   0.3868     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_species_class)
   R2_marginal R2_conditional 
    0.02688044     0.36823812 
Code
orchard_plot(mod_sex_species_class, mod = "sex_species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 0)

Figure S21. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_species_class <- all_models(mod_sex_species_class,  mod = "sex_species_class", type = "normal")

# saving tab_sex_species_class as RDS
saveRDS(tab_sex_species_class, here("Rdata", "tab_sex_species_class.RDS"))

Table S14. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_species_class <-readRDS(here("Rdata", "tab_sex_species_class.RDS"))

tab_sex_species_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_Aves -0.001 -0.116 0.114 0.987 -0.777 0.775
B_Insecta -0.146 -0.689 0.397 0.597 -1.086 0.794
B_Mammalia -0.180 -0.312 -0.047 0.008 -0.959 0.599
B_Reptilia 0.158 -0.381 0.696 0.565 -0.780 1.095
F_Actinopterygii -0.069 -0.477 0.339 0.741 -0.938 0.800
F_Arachnida 0.040 -0.514 0.594 0.887 -0.906 0.987
F_Aves -0.036 -0.124 0.051 0.417 -0.809 0.736
F_Insecta 0.240 -0.128 0.607 0.201 -0.611 1.091
F_Mammalia -0.019 -0.124 0.087 0.731 -0.793 0.756
F_Reptilia 0.071 -0.238 0.380 0.651 -0.756 0.898
M_Actinopterygii 0.068 -0.339 0.475 0.742 -0.800 0.937
M_Aves -0.035 -0.123 0.054 0.441 -0.807 0.738
M_Insecta -0.143 -0.888 0.602 0.707 -1.212 0.927
M_Mammalia -0.016 -0.129 0.097 0.778 -0.792 0.760
M_Reptilia -0.024 -0.481 0.433 0.918 -0.917 0.869
B_Aves-B_Insecta -0.145 -0.700 0.409 0.607 -1.092 0.801
B_Aves-B_Mammalia -0.179 -0.346 -0.012 0.036 -0.964 0.607
B_Aves-B_Reptilia 0.159 -0.388 0.706 0.569 -0.784 1.101
B_Aves-F_Actinopterygii -0.068 -0.489 0.353 0.752 -0.943 0.807
B_Aves-F_Arachnida 0.041 -0.523 0.605 0.886 -0.912 0.994
B_Aves-F_Aves -0.035 -0.165 0.094 0.594 -0.814 0.743
B_Aves-F_Insecta 0.241 -0.141 0.623 0.217 -0.617 1.098
B_Aves-F_Mammalia -0.018 -0.165 0.130 0.815 -0.799 0.764
B_Aves-F_Reptilia 0.072 -0.253 0.396 0.663 -0.761 0.905
B_Aves-M_Actinopterygii 0.069 -0.351 0.490 0.746 -0.806 0.944
B_Aves-M_Aves -0.034 -0.163 0.096 0.611 -0.812 0.745
B_Aves-M_Insecta -0.142 -0.894 0.611 0.712 -1.217 0.933
B_Aves-M_Mammalia -0.015 -0.169 0.138 0.845 -0.798 0.767
B_Aves-M_Reptilia -0.023 -0.491 0.445 0.923 -0.922 0.876
B_Insecta-B_Mammalia -0.033 -0.591 0.524 0.906 -0.982 0.915
B_Insecta-B_Reptilia 0.304 -0.459 1.067 0.435 -0.779 1.386
B_Insecta-F_Actinopterygii 0.077 -0.601 0.756 0.823 -0.947 1.102
B_Insecta-F_Arachnida 0.186 -0.589 0.962 0.637 -0.904 1.277
B_Insecta-F_Aves 0.110 -0.439 0.659 0.694 -0.834 1.054
B_Insecta-F_Insecta 0.386 -0.269 1.041 0.248 -0.623 1.395
B_Insecta-F_Mammalia 0.128 -0.424 0.680 0.650 -0.818 1.073
B_Insecta-F_Reptilia 0.217 -0.406 0.841 0.494 -0.771 1.206
B_Insecta-M_Actinopterygii 0.215 -0.463 0.892 0.534 -0.809 1.239
B_Insecta-M_Aves 0.112 -0.438 0.661 0.690 -0.832 1.056
B_Insecta-M_Insecta 0.004 -0.918 0.925 0.994 -1.196 1.203
B_Insecta-M_Mammalia 0.130 -0.423 0.684 0.645 -0.816 1.076
B_Insecta-M_Reptilia 0.122 -0.587 0.831 0.735 -0.922 1.167
B_Mammalia-B_Reptilia 0.337 -0.210 0.885 0.227 -0.605 1.280
B_Mammalia-F_Actinopterygii 0.111 -0.312 0.534 0.607 -0.765 0.987
B_Mammalia-F_Arachnida 0.220 -0.347 0.786 0.446 -0.734 1.174
B_Mammalia-F_Aves 0.144 -0.004 0.291 0.056 -0.638 0.925
B_Mammalia-F_Insecta 0.419 0.035 0.804 0.033 -0.439 1.278
B_Mammalia-F_Mammalia 0.161 0.020 0.303 0.026 -0.619 0.942
B_Mammalia-F_Reptilia 0.251 -0.075 0.577 0.132 -0.583 1.085
B_Mammalia-M_Actinopterygii 0.248 -0.175 0.671 0.250 -0.628 1.124
B_Mammalia-M_Aves 0.145 -0.002 0.293 0.054 -0.636 0.927
B_Mammalia-M_Insecta 0.037 -0.717 0.791 0.923 -1.039 1.113
B_Mammalia-M_Mammalia 0.163 0.015 0.312 0.031 -0.618 0.945
B_Mammalia-M_Reptilia 0.156 -0.314 0.625 0.515 -0.744 1.055
B_Reptilia-F_Actinopterygii -0.226 -0.897 0.444 0.507 -1.245 0.793
B_Reptilia-F_Arachnida -0.117 -0.887 0.652 0.764 -1.204 0.969
B_Reptilia-F_Aves -0.194 -0.735 0.347 0.482 -1.133 0.745
B_Reptilia-F_Insecta 0.082 -0.565 0.729 0.803 -0.922 1.086
B_Reptilia-F_Mammalia -0.176 -0.718 0.366 0.524 -1.116 0.764
B_Reptilia-F_Reptilia -0.087 -0.653 0.479 0.764 -1.040 0.867
B_Reptilia-M_Actinopterygii -0.089 -0.760 0.581 0.794 -1.109 0.930
B_Reptilia-M_Aves -0.192 -0.733 0.349 0.486 -1.131 0.747
B_Reptilia-M_Insecta -0.300 -1.217 0.616 0.520 -1.496 0.895
B_Reptilia-M_Mammalia -0.174 -0.718 0.370 0.530 -1.115 0.767
B_Reptilia-M_Reptilia -0.182 -0.860 0.497 0.599 -1.206 0.843
F_Actinopterygii-F_Arachnida 0.109 -0.577 0.795 0.755 -0.920 1.138
F_Actinopterygii-F_Aves 0.033 -0.381 0.446 0.877 -0.839 0.904
F_Actinopterygii-F_Insecta 0.308 -0.237 0.854 0.267 -0.633 1.250
F_Actinopterygii-F_Mammalia 0.050 -0.366 0.466 0.812 -0.823 0.923
F_Actinopterygii-F_Reptilia 0.140 -0.366 0.646 0.588 -0.779 1.059
F_Actinopterygii-M_Actinopterygii 0.137 -0.324 0.598 0.559 -0.758 1.032
F_Actinopterygii-M_Aves 0.034 -0.379 0.448 0.871 -0.838 0.906
F_Actinopterygii-M_Insecta -0.074 -0.921 0.773 0.864 -1.217 1.069
F_Actinopterygii-M_Mammalia 0.053 -0.366 0.471 0.805 -0.821 0.927
F_Actinopterygii-M_Reptilia 0.045 -0.564 0.653 0.885 -0.935 1.024
F_Arachnida-F_Aves -0.076 -0.635 0.483 0.789 -1.026 0.873
F_Arachnida-F_Insecta 0.200 -0.463 0.862 0.555 -0.814 1.213
F_Arachnida-F_Mammalia -0.059 -0.620 0.502 0.837 -1.009 0.892
F_Arachnida-F_Reptilia 0.031 -0.600 0.662 0.924 -0.963 1.024
F_Arachnida-M_Actinopterygii 0.028 -0.657 0.714 0.936 -1.001 1.057
F_Arachnida-M_Aves -0.075 -0.634 0.484 0.793 -1.024 0.875
F_Arachnida-M_Insecta -0.183 -1.110 0.744 0.699 -1.386 1.021
F_Arachnida-M_Mammalia -0.056 -0.619 0.506 0.844 -1.008 0.895
F_Arachnida-M_Reptilia -0.064 -0.780 0.651 0.860 -1.114 0.985
F_Aves-F_Insecta 0.276 -0.098 0.650 0.148 -0.578 1.130
F_Aves-F_Mammalia 0.018 -0.108 0.143 0.782 -0.760 0.795
F_Aves-F_Reptilia 0.107 -0.208 0.422 0.504 -0.722 0.937
F_Aves-M_Actinopterygii 0.104 -0.309 0.518 0.620 -0.767 0.976
F_Aves-M_Aves 0.002 -0.074 0.077 0.967 -0.770 0.773
F_Aves-M_Insecta -0.107 -0.855 0.642 0.780 -1.179 0.965
F_Aves-M_Mammalia 0.020 -0.112 0.152 0.768 -0.759 0.799
F_Aves-M_Reptilia 0.012 -0.449 0.474 0.959 -0.884 0.908
F_Insecta-F_Mammalia -0.258 -0.635 0.118 0.179 -1.113 0.597
F_Insecta-F_Reptilia -0.169 -0.643 0.305 0.485 -1.071 0.733
F_Insecta-M_Actinopterygii -0.171 -0.717 0.374 0.537 -1.113 0.770
F_Insecta-M_Aves -0.274 -0.648 0.100 0.151 -1.128 0.580
F_Insecta-M_Insecta -0.382 -1.165 0.400 0.338 -1.479 0.714
F_Insecta-M_Mammalia -0.256 -0.635 0.123 0.185 -1.112 0.600
F_Insecta-M_Reptilia -0.264 -0.846 0.318 0.374 -1.227 0.700
F_Mammalia-F_Reptilia 0.090 -0.228 0.407 0.580 -0.741 0.920
F_Mammalia-M_Actinopterygii 0.087 -0.329 0.503 0.682 -0.786 0.960
F_Mammalia-M_Aves -0.016 -0.141 0.109 0.802 -0.794 0.762
F_Mammalia-M_Insecta -0.124 -0.874 0.626 0.745 -1.197 0.949
F_Mammalia-M_Mammalia 0.002 -0.099 0.103 0.965 -0.772 0.776
F_Mammalia-M_Reptilia -0.006 -0.469 0.458 0.981 -0.902 0.891
F_Reptilia-M_Actinopterygii -0.003 -0.509 0.503 0.992 -0.922 0.917
F_Reptilia-M_Aves -0.106 -0.420 0.209 0.511 -0.935 0.724
F_Reptilia-M_Insecta -0.214 -1.017 0.590 0.602 -1.325 0.898
F_Reptilia-M_Mammalia -0.087 -0.407 0.233 0.593 -0.919 0.744
F_Reptilia-M_Reptilia -0.095 -0.542 0.352 0.676 -0.983 0.793
M_Actinopterygii-M_Aves -0.103 -0.516 0.311 0.625 -0.975 0.769
M_Actinopterygii-M_Insecta -0.211 -1.058 0.636 0.625 -1.354 0.932
M_Actinopterygii-M_Mammalia -0.085 -0.503 0.334 0.691 -0.959 0.789
M_Actinopterygii-M_Reptilia -0.092 -0.701 0.516 0.766 -1.072 0.887
M_Aves-M_Insecta -0.108 -0.857 0.640 0.777 -1.180 0.964
M_Aves-M_Mammalia 0.018 -0.114 0.151 0.786 -0.761 0.797
M_Aves-M_Reptilia 0.010 -0.451 0.472 0.965 -0.885 0.906
M_Insecta-M_Mammalia 0.126 -0.625 0.878 0.741 -0.948 1.201
M_Insecta-M_Reptilia 0.119 -0.753 0.990 0.789 -1.043 1.280
M_Mammalia-M_Reptilia -0.008 -0.473 0.457 0.974 -0.905 0.890
Code
dat$sex_whose_fitness <- as.factor(paste(dat$sex, dat$whose_fitness ,sep = "_"))

mod_sex_whose_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_whose_fitness - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_whose_fitness)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.8293   481.6587   499.6587   540.7587   499.9155   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0726  0.2695    717     no         effectID 
sigma^2.2  0.0243  0.1560    210     no          paperID 
sigma^2.3  0.0135  0.1164    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 711) = 15581.8457, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 711) = 1.2488, p-val = 0.2792

Model Results:

                                    estimate      se     tval   df    pval 
sex_whose_fitnessBoth_descendant      0.0376  0.1083   0.3470  711  0.7287 
sex_whose_fitnessBoth_focal          -0.0879  0.0422  -2.0833  711  0.0376 
sex_whose_fitnessFemale_descendant   -0.0115  0.0491  -0.2352  711  0.8141 
sex_whose_fitnessFemale_focal        -0.0027  0.0340  -0.0789  711  0.9372 
sex_whose_fitnessMale_descendant      0.0592  0.0553   1.0702  711  0.2849 
sex_whose_fitnessMale_focal          -0.0120  0.0350  -0.3444  711  0.7307 
                                      ci.lb    ci.ub    
sex_whose_fitnessBoth_descendant    -0.1751   0.2503    
sex_whose_fitnessBoth_focal         -0.1707  -0.0051  * 
sex_whose_fitnessFemale_descendant  -0.1078   0.0848    
sex_whose_fitnessFemale_focal       -0.0694   0.0640    
sex_whose_fitnessMale_descendant    -0.0494   0.1677    
sex_whose_fitnessMale_focal         -0.0807   0.0566    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_whose_fitness)
   R2_marginal R2_conditional 
    0.01133881     0.35022215 
Code
orchard_plot(mod_sex_whose_fitness, mod = "sex_whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 0)

Figure S22. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_whose_fitness <- all_models(mod_sex_whose_fitness,  mod = "sex_whose_fitness", type = "normal")

# saving tab_sex_whose_fitness as RDS
saveRDS(tab_sex_whose_fitness, here("Rdata", "tab_sex_whose_fitness.RDS"))

Table S15. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_whose_fitness <-readRDS(here("Rdata", "tab_sex_whose_fitness.RDS"))

tab_sex_whose_fitness
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_individual -0.078 -0.168 0.011 0.086 -0.849 0.693
B_offspring 0.020 -0.276 0.316 0.894 -0.801 0.841
F_individual -0.011 -0.084 0.062 0.760 -0.781 0.758
F_offspring -0.021 -0.129 0.087 0.704 -0.794 0.753
M_individual -0.034 -0.109 0.041 0.376 -0.803 0.736
M_offspring 0.061 -0.063 0.185 0.333 -0.715 0.837
B_individual-B_offspring 0.098 -0.198 0.395 0.515 -0.723 0.920
B_individual-F_individual 0.067 -0.028 0.162 0.168 -0.705 0.839
B_individual-F_offspring 0.057 -0.067 0.182 0.365 -0.719 0.833
B_individual-M_individual 0.045 -0.053 0.142 0.370 -0.728 0.817
B_individual-M_offspring 0.140 0.002 0.277 0.047 -0.639 0.918
B_offspring-F_individual -0.031 -0.329 0.266 0.836 -0.853 0.790
B_offspring-F_offspring -0.041 -0.349 0.268 0.795 -0.867 0.785
B_offspring-M_individual -0.054 -0.352 0.244 0.723 -0.876 0.768
B_offspring-M_offspring 0.041 -0.273 0.355 0.797 -0.787 0.869
F_individual-F_offspring -0.010 -0.112 0.093 0.855 -0.782 0.763
F_individual-M_individual -0.022 -0.087 0.043 0.499 -0.791 0.746
F_individual-M_offspring 0.073 -0.046 0.192 0.232 -0.703 0.848
F_offspring-M_individual -0.013 -0.118 0.092 0.810 -0.786 0.760
F_offspring-M_offspring 0.082 -0.051 0.215 0.227 -0.695 0.860
M_individual-M_offspring 0.095 -0.025 0.215 0.121 -0.680 0.870
Code
dat$sex_fitness_higher_level <- as.factor(paste(dat$sex, dat$fitness_higher_level ,sep = "_"))

mod_sex_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_fitness_higher_level - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_fitness)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.1219   480.2438   502.2438   552.4462   502.6226   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0722  0.2687    717     no         effectID 
sigma^2.2  0.0265  0.1627    210     no          paperID 
sigma^2.3  0.0115  0.1072    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 709) = 15509.7239, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 709) = 1.2527, p-val = 0.2655

Model Results:

                                                 estimate      se     tval   df 
sex_fitness_higher_levelBoth_reproduction         -0.0199  0.0622  -0.3209  709 
sex_fitness_higher_levelBoth_survival             -0.1022  0.0469  -2.1812  709 
sex_fitness_higher_levelFemale_indirect fitness   -0.2290  0.2110  -1.0852  709 
sex_fitness_higher_levelFemale_reproduction        0.0091  0.0358   0.2543  709 
sex_fitness_higher_levelFemale_survival           -0.0267  0.0390  -0.6863  709 
sex_fitness_higher_levelMale_indirect fitness     -0.2261  0.1728  -1.3087  709 
sex_fitness_higher_levelMale_reproduction         -0.0135  0.0386  -0.3486  709 
sex_fitness_higher_levelMale_survival              0.0185  0.0404   0.4568  709 
                                                   pval    ci.lb    ci.ub    
sex_fitness_higher_levelBoth_reproduction        0.7484  -0.1420   0.1021    
sex_fitness_higher_levelBoth_survival            0.0295  -0.1942  -0.0102  * 
sex_fitness_higher_levelFemale_indirect fitness  0.2782  -0.6433   0.1853    
sex_fitness_higher_levelFemale_reproduction      0.7994  -0.0611   0.0793    
sex_fitness_higher_levelFemale_survival          0.4927  -0.1032   0.0498    
sex_fitness_higher_levelMale_indirect fitness    0.1911  -0.5654   0.1131    
sex_fitness_higher_levelMale_reproduction        0.7275  -0.0893   0.0623    
sex_fitness_higher_levelMale_survival            0.6479  -0.0609   0.0978    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_fitness)
   R2_marginal R2_conditional 
    0.01409984     0.35376375 
Code
orchard_plot(mod_sex_fitness, mod = "sex_fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 0)

Figure S23. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_fitness <- all_models(mod_sex_fitness,  mod = "sex_fitness_higher_level", type = "normal")

# saving tab_sex_fitness as RDS
saveRDS(tab_sex_fitness, here("Rdata", "tab_sex_fitness.RDS"))

Table S16. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_fitness <-readRDS(here("Rdata", "tab_sex_fitness.RDS"))

tab_sex_fitness
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_reproduction -0.023 -0.164 0.118 0.749 -0.800 0.754
B_survival -0.095 -0.196 0.006 0.066 -0.865 0.676
F_indirect fitness -0.231 -0.732 0.270 0.365 -1.144 0.682
F_reproduction -0.004 -0.081 0.073 0.922 -0.771 0.764
F_survival -0.026 -0.112 0.061 0.558 -0.794 0.743
M_indirect fitness -0.230 -0.638 0.178 0.269 -1.096 0.636
M_reproduction -0.028 -0.112 0.056 0.513 -0.796 0.740
M_survival -0.001 -0.089 0.088 0.988 -0.769 0.768
B_reproduction-B_survival -0.072 -0.231 0.088 0.378 -0.852 0.708
B_reproduction-F_indirect fitness -0.208 -0.726 0.310 0.431 -1.131 0.715
B_reproduction-F_reproduction 0.019 -0.132 0.170 0.803 -0.759 0.798
B_reproduction-F_survival -0.003 -0.159 0.153 0.972 -0.782 0.777
B_reproduction-M_indirect fitness -0.207 -0.636 0.223 0.345 -1.083 0.669
B_reproduction-M_reproduction -0.005 -0.159 0.149 0.949 -0.784 0.774
B_reproduction-M_survival 0.022 -0.135 0.180 0.781 -0.757 0.802
B_survival-F_indirect fitness -0.136 -0.644 0.371 0.598 -1.053 0.781
B_survival-F_reproduction 0.091 -0.017 0.199 0.098 -0.680 0.862
B_survival-F_survival 0.069 -0.045 0.183 0.233 -0.703 0.841
B_survival-M_indirect fitness -0.135 -0.552 0.282 0.525 -1.005 0.735
B_survival-M_reproduction 0.067 -0.046 0.180 0.247 -0.705 0.839
B_survival-M_survival 0.094 -0.023 0.211 0.114 -0.678 0.867
F_indirect fitness-F_reproduction 0.227 -0.275 0.729 0.375 -0.687 1.141
F_indirect fitness-F_survival 0.205 -0.297 0.708 0.423 -0.709 1.119
F_indirect fitness-M_indirect fitness 0.001 -0.602 0.604 0.997 -0.972 0.974
F_indirect fitness-M_reproduction 0.203 -0.300 0.706 0.428 -0.711 1.118
F_indirect fitness-M_survival 0.230 -0.272 0.733 0.369 -0.684 1.145
F_reproduction-F_survival -0.022 -0.104 0.060 0.600 -0.790 0.746
F_reproduction-M_indirect fitness -0.226 -0.636 0.184 0.279 -1.093 0.641
F_reproduction-M_reproduction -0.024 -0.102 0.054 0.542 -0.792 0.743
F_reproduction-M_survival 0.003 -0.085 0.091 0.944 -0.766 0.772
F_survival-M_indirect fitness -0.204 -0.615 0.207 0.330 -1.071 0.663
F_survival-M_reproduction -0.002 -0.093 0.089 0.962 -0.771 0.767
F_survival-M_survival 0.025 -0.064 0.114 0.579 -0.744 0.794
M_indirect fitness-M_reproduction 0.202 -0.208 0.612 0.334 -0.665 1.069
M_indirect fitness-M_survival 0.229 -0.182 0.640 0.274 -0.638 1.096
M_reproduction-M_survival 0.027 -0.067 0.122 0.570 -0.742 0.797
Code
dat$sex_dispersal_type <- as.factor(paste(dat$sex, dat$dispersal_type ,sep = "_"))

mod_sex_dispersal_type <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_dispersal_type - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_dispersal_type)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.8846   481.7693   505.7693   560.5186   506.2182   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0728  0.2698    717     no         effectID 
sigma^2.2  0.0320  0.1788    210     no          paperID 
sigma^2.3  0.0059  0.0770    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 708) = 15380.1268, p-val < .0001

Test of Moderators (coefficients 1:9):
F(df1 = 9, df2 = 708) = 0.7656, p-val = 0.6485

Model Results:

                                   estimate      se     tval   df    pval 
sex_dispersal_typeBoth_Both         -0.1156  0.1200  -0.9635  708  0.3356 
sex_dispersal_typeBoth_breeding      0.0162  0.0829   0.1959  708  0.8448 
sex_dispersal_typeBoth_natal        -0.0949  0.0464  -2.0444  708  0.0413 
sex_dispersal_typeFemale_Both       -0.0877  0.1205  -0.7281  708  0.4668 
sex_dispersal_typeFemale_breeding   -0.0069  0.0489  -0.1409  708  0.8880 
sex_dispersal_typeFemale_natal       0.0014  0.0354   0.0397  708  0.9683 
sex_dispersal_typeMale_Both         -0.0066  0.1371  -0.0479  708  0.9618 
sex_dispersal_typeMale_breeding     -0.0257  0.0524  -0.4912  708  0.6234 
sex_dispersal_typeMale_natal         0.0127  0.0373   0.3414  708  0.7329 
                                     ci.lb    ci.ub    
sex_dispersal_typeBoth_Both        -0.3511   0.1199    
sex_dispersal_typeBoth_breeding    -0.1465   0.1790    
sex_dispersal_typeBoth_natal       -0.1860  -0.0038  * 
sex_dispersal_typeFemale_Both      -0.3243   0.1488    
sex_dispersal_typeFemale_breeding  -0.1030   0.0892    
sex_dispersal_typeFemale_natal     -0.0681   0.0709    
sex_dispersal_typeMale_Both        -0.2757   0.2625    
sex_dispersal_typeMale_breeding    -0.1286   0.0771    
sex_dispersal_typeMale_natal       -0.0604   0.0859    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_dispersal_type)
   R2_marginal R2_conditional 
    0.01210331     0.35021950 
Code
orchard_plot(mod_sex_dispersal_type, mod = "sex_dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, trunk.size = 0.5, angle = 0)

Figure S24. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types combined). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_dispersal_type <- all_models(mod_sex_dispersal_type,  mod = "sex_dispersal_type", type = "normal")

# saving tab_sex_dispersal_type as RDS
saveRDS(tab_sex_dispersal_type, here("Rdata", "tab_sex_dispersal_type.RDS"))

Table S17. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types grouped) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_dispersal_type <-readRDS(here("Rdata", "tab_sex_dispersal_type.RDS"))

tab_sex_dispersal_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_B -0.096 -0.321 0.129 0.403 -0.892 0.699
B_breeding -0.004 -0.200 0.192 0.967 -0.792 0.784
B_natal -0.083 -0.184 0.017 0.103 -0.853 0.686
F_B -0.081 -0.301 0.138 0.468 -0.875 0.713
F_breeding -0.004 -0.110 0.102 0.940 -0.774 0.766
F_natal -0.008 -0.086 0.070 0.840 -0.775 0.759
M_B -0.108 -0.350 0.133 0.379 -0.909 0.692
M_breeding -0.022 -0.137 0.094 0.712 -0.793 0.750
M_natal -0.004 -0.086 0.077 0.922 -0.771 0.763
B_B-B_breeding 0.092 -0.203 0.387 0.540 -0.726 0.910
B_B-B_natal 0.013 -0.223 0.249 0.916 -0.786 0.811
B_B-F_B 0.015 -0.288 0.317 0.923 -0.806 0.836
B_B-F_breeding 0.092 -0.150 0.334 0.456 -0.708 0.892
B_B-F_natal 0.088 -0.142 0.318 0.453 -0.709 0.885
B_B-M_B -0.012 -0.330 0.305 0.939 -0.839 0.814
B_B-M_breeding 0.074 -0.172 0.321 0.554 -0.727 0.876
B_B-M_natal 0.092 -0.139 0.323 0.435 -0.705 0.889
B_breeding-B_natal -0.079 -0.292 0.133 0.464 -0.871 0.713
B_breeding-F_B -0.077 -0.368 0.214 0.603 -0.893 0.739
B_breeding-F_breeding 0.000 -0.214 0.214 1.000 -0.792 0.792
B_breeding-F_natal -0.004 -0.208 0.200 0.970 -0.794 0.786
B_breeding-M_B -0.104 -0.413 0.204 0.506 -0.927 0.718
B_breeding-M_breeding -0.018 -0.238 0.202 0.875 -0.812 0.776
B_breeding-M_natal 0.000 -0.206 0.206 1.000 -0.790 0.790
B_natal-F_B 0.002 -0.232 0.237 0.985 -0.796 0.800
B_natal-F_breeding 0.079 -0.051 0.209 0.232 -0.695 0.853
B_natal-F_natal 0.075 -0.032 0.183 0.169 -0.695 0.846
B_natal-M_B -0.025 -0.281 0.231 0.848 -0.830 0.780
B_natal-M_breeding 0.062 -0.078 0.202 0.386 -0.714 0.837
B_natal-M_natal 0.079 -0.031 0.190 0.158 -0.692 0.850
F_B-F_breeding 0.077 -0.160 0.314 0.523 -0.722 0.876
F_B-F_natal 0.073 -0.153 0.300 0.526 -0.723 0.869
F_B-M_B -0.027 -0.288 0.234 0.837 -0.834 0.779
F_B-M_breeding 0.059 -0.183 0.302 0.630 -0.741 0.860
F_B-M_natal 0.077 -0.150 0.304 0.506 -0.719 0.873
F_breeding-F_natal -0.004 -0.113 0.105 0.944 -0.775 0.767
F_breeding-M_B -0.104 -0.362 0.154 0.427 -0.910 0.701
F_breeding-M_breeding -0.018 -0.130 0.095 0.758 -0.789 0.754
F_breeding-M_natal 0.000 -0.112 0.112 1.000 -0.771 0.771
F_natal-M_B -0.100 -0.349 0.148 0.427 -0.903 0.702
F_natal-M_breeding -0.014 -0.134 0.107 0.823 -0.786 0.759
F_natal-M_natal 0.004 -0.068 0.076 0.916 -0.762 0.770
M_B-M_breeding 0.087 -0.176 0.350 0.517 -0.720 0.894
M_B-M_natal 0.104 -0.145 0.354 0.412 -0.698 0.907
M_breeding-M_natal 0.018 -0.106 0.142 0.780 -0.755 0.791

5 Species-level models

There are two steps, which can be done using phylo_model().

5.1 Estimate species-level effect sizes.

Code
# function to estimate species-level effect sizes
phylo_model <- function(data = dat, tree, correlation = 0.5) {
  
  # accounting for within-species correlation
  Vm <- vcalc(vi = data$vi, cluster = data$shared_group, obs = data$effectID, rho = correlation)
  mod <- rma.mv(yi = yi, 
              V = Vm,
              mod = ~ 1,
              data = data,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | phylogeny),
              #R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)
  
  
  # get random effects
  # dat3 <- dat2[c("phylogeny", "yi", "vi")]
  REs <- ranef(mod)$phylogeny
  
  # calculate means and CIs
  dat <- data.frame(phylogeny = rownames(REs),
                    estimate = mod$beta[1] + REs$intrcpt,
                    se = sqrt((mod$se[1])^2 + (REs$se)^2)) %>%
    mutate(lb = estimate - se * qnorm(1 - 0.05/2), 
           ub = estimate + se * qnorm(1 - 0.05/2)) %>%
    arrange(phylogeny)
  # get species class data
  species_class <- dplyr::distinct(data, phylogeny, .keep_all = TRUE) %>%
    dplyr::select(species_class, phylogeny)
  
  # extract tip labels from the tree and merge with species class
  tip.label <- data.frame(phylogeny = tree$tip.label)  # extract tip label
  species_class2 <- left_join(tip.label, species_class, by = "phylogeny")
  
  return(list(data = dat, class = species_class2))  # return processed data and species class
}

# estimate species-level effect sizes
# you can specify your own correlation coefficient via `correlation`. For illustration purpose, I used 0.5.
result <- phylo_model(data = dat, tree = tree, correlation = 0.5)

# find the significant species
result$data <- result$data %>%
  mutate(
    p_value = 2 * (1 - pnorm(abs(estimate / se))),
    significance = ifelse(p_value < 0.05, "Significant", "NonSignificant")
  )
result$data <- result$data %>%
  mutate(significance_label = ifelse(p_value < 0.05, "*", NA))

5.2 Plot tree and add species-level effect size estimates

Code
# plot tree
p1 <- ggtree(tree, layout = "rectangular", cex = 0.3)

p2 <- p1 %<+% result$class + geom_tiplab(aes(color = species_class), size = 3, align = T, offset = 0.05) + geom_tippoint(aes(color = species_class)) + guides(color = "none")

p3 <- p2 + 
  geom_facet(
    panel = "Effect size", 
    data = result$data, 
    geom = ggstance::geom_pointrangeh,
    mapping = aes(x = estimate, xmin = lb, xmax = ub, color = species_class)
  ) +
  geom_facet(
    panel = "Effect size",
    data = result$data,
    geom = geom_text,
    mapping = aes(x = lb, label = significance_label), 
    hjust = 1.5, size = 3, color = "red"              # adjust position of star *
  ) +
  geom_facet(
    panel = "Effect size",
    data = result$data, 
    geom = geom_vline,
    mapping = aes(xintercept = 0),                    
    linetype = "dashed", color = "gray"
  ) +
  theme_tree2()
Code
# adjust widths of the plot facets as necessary to improve visualization 

facet_widths(p3, c(Tree = 0.6, `Effect size` = 0.4))

Figure S25. Phylogenetic tree (left panel) and species-level effect size estimates (right panel). Significant relationships are indicated with a red asterisk. Dispersal was positive for ruffed grouse (Bonasa umbellus; Zr = 0.42, 95% CIs: 0.16 to 0.68) and negative for brown jays (Cyanocorax morio; Zr = -0.30, 95% CIs: -0.58 to -0.01), Tasmanian nativehens (Gallinula mortierii; Zr = -0.34, 95% CIs: -0.60 to X), and cougars (Puma concolor; Zr = -0.61, 95% CIs: -0.90 to -0.32).

6 Multi-moderator-regression model

Code
# used "maximum likelihood" for model selection

#######################
# Multi-variable models
#######################

mod_full <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1 + sex +
                age_class_clean + whose_fitness + 
                fitness_higher_level +
                dispersal_type + dispersal_phase,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              method = "ML",
              test = "t",
              sparse = TRUE)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 717; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-237.7901  2915.8747   503.5803   567.6313   504.1785   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0714  0.2672    717     no         effectID 
sigma^2.2  0.0198  0.1407    210     no          paperID 
sigma^2.3  0.0180  0.1343    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 706) = 15385.8283, p-val < .0001

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 706) = 1.2564, p-val = 0.2515

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.2536  0.1723  -1.4719  706  0.1415 
sexFemale                           0.0583  0.0433   1.3467  706  0.1785 
sexMale                             0.0662  0.0438   1.5114  706  0.1311 
age_class_cleanjuvenile             0.0830  0.0551   1.5058  706  0.1326 
age_class_cleanmix                 -0.0569  0.0604  -0.9423  706  0.3464 
whose_fitnessfocal                 -0.0534  0.0406  -1.3144  706  0.1892 
fitness_higher_levelreproduction    0.2270  0.1429   1.5887  706  0.1126 
fitness_higher_levelsurvival        0.1876  0.1448   1.2960  706  0.1954 
dispersal_typebreeding             -0.0059  0.0885  -0.0672  706  0.9464 
dispersal_typenatal                -0.0053  0.0810  -0.0656  706  0.9477 
dispersal_phasesettlement           0.0311  0.0491   0.6336  706  0.5265 
                                    ci.lb   ci.ub    
intrcpt                           -0.5920  0.0847    
sexFemale                         -0.0267  0.1432    
sexMale                           -0.0198  0.1522    
age_class_cleanjuvenile           -0.0252  0.1913    
age_class_cleanmix                -0.1754  0.0616    
whose_fitnessfocal                -0.1332  0.0264    
fitness_higher_levelreproduction  -0.0535  0.5075    
fitness_higher_levelsurvival      -0.0966  0.4719    
dispersal_typebreeding            -0.1796  0.1677    
dispersal_typenatal               -0.1644  0.1538    
dispersal_phasesettlement         -0.0653  0.1276    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_full)*100, 2)
   R2_marginal R2_conditional 
          2.24          36.09 
Code
# multi-model selection
candidates <- dredge(mod_full, trace = 2)

# saving candidates as RDS
saveRDS(candidates, here("Rdata", "candidates.RDS"))
Code
candidates <-readRDS(here("Rdata", "candidates.RDS"))

candidates
Global model call: rma.mv(yi = yi, V = VCV, mods = ~1 + sex + age_class_clean + 
    whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, 
    random = list(~1 | effectID, ~1 | paperID, ~1 | species_cleaned), 
    data = dat, method = "ML", test = "t", sparse = TRUE)
---
Model selection table 
   (Int) age_cls_cln dsp_phs dsp_typ ftn_hgh_lvl sex whs_ftn df   logLik  AICc
2      +           +                                          6 -276.828 565.8
4      +           +       +                                  7 -275.933 566.0
34     +           +                                       +  7 -276.278 566.7
36     +           +       +                               +  8 -275.479 567.2
1      +                                                      4 -279.768 567.6
42     +           +                           +           +  9 -274.701 567.7
10     +           +                           +              8 -275.782 567.8
18     +           +                               +          8 -275.927 568.1
12     +           +       +                   +              9 -275.053 568.4
33     +                                                   +  5 -279.217 568.5
3      +                   +                                  5 -279.276 568.6
44     +           +       +                   +           + 10 -274.319 569.0
20     +           +       +                       +          9 -275.359 569.0
50     +           +                               +       +  9 -275.518 569.3
6      +           +               +                          8 -276.641 569.5
9      +                                       +              6 -278.813 569.8
35     +                   +                               +  6 -278.815 569.8
41     +                                       +           +  7 -277.874 569.9
8      +           +       +       +                          9 -275.830 569.9
17     +                                           +          6 -278.913 570.0
26     +           +                           +   +         10 -274.949 570.2
52     +           +       +                       +       + 10 -275.002 570.3
38     +           +               +                       +  9 -276.064 570.4
5      +                           +                          6 -279.306 570.7
58     +           +                           +   +       + 11 -274.187 570.8
40     +           +       +       +                       + 10 -275.349 571.0
11     +                   +                   +              7 -278.449 571.1
49     +                                           +       +  7 -278.515 571.2
28     +           +       +                   +   +         11 -274.474 571.3
46     +           +               +           +           + 11 -274.557 571.5
19     +                   +                       +          7 -278.735 571.6
14     +           +               +           +             10 -275.674 571.7
37     +                           +                       +  7 -278.774 571.7
43     +                   +                   +           +  8 -277.769 571.8
22     +           +               +               +         10 -275.770 571.9
7      +                   +       +                          7 -278.938 572.0
25     +                                       +   +          8 -278.038 572.3
60     +           +       +                   +   +       + 12 -273.929 572.3
16     +           +       +       +           +             11 -274.996 572.4
48     +           +       +       +           +           + 12 -274.220 572.9
24     +           +       +       +               +         11 -275.261 572.9
51     +                   +                       +       +  8 -278.367 572.9
54     +           +               +               +       + 11 -275.338 573.1
57     +                                       +   +       +  9 -277.439 573.1
13     +                           +           +              8 -278.487 573.2
39     +                   +       +                       +  8 -278.490 573.2
21     +                           +               +          8 -278.571 573.4
45     +                           +           +           +  9 -277.602 573.5
27     +                   +                   +   +          9 -277.898 574.1
30     +           +               +           +   +         12 -274.856 574.2
56     +           +       +       +               +       + 12 -274.881 574.2
53     +                           +               +       +  9 -278.186 574.6
15     +                   +       +           +              9 -278.197 574.7
62     +           +               +           +   +       + 13 -274.062 574.7
59     +                   +                   +   +       + 10 -277.406 575.1
23     +                   +       +               +          9 -278.444 575.2
32     +           +       +       +           +   +         13 -274.418 575.4
47     +                   +       +           +           + 10 -277.534 575.4
29     +                           +           +   +         10 -277.795 575.9
64     +           +       +       +           +   +       + 14 -273.837 576.3
55     +                   +       +               +       + 10 -278.088 576.5
61     +                           +           +   +       + 11 -277.225 576.8
31     +                   +       +           +   +         11 -277.687 577.8
63     +                   +       +           +   +       + 12 -277.207 578.9
   delta weight
2   0.00  0.122
4   0.25  0.107
34  0.94  0.076
36  1.39  0.061
1   1.81  0.049
42  1.89  0.047
10  2.00  0.045
18  2.29  0.039
12  2.59  0.033
33  2.74  0.031
3   2.86  0.029
44  3.19  0.025
20  3.21  0.024
50  3.52  0.021
6   3.72  0.019
9   3.97  0.017
35  3.97  0.017
41  4.13  0.015
8   4.15  0.015
17  4.17  0.015
26  4.45  0.013
52  4.55  0.012
38  4.62  0.012
5   4.96  0.010
58  4.99  0.010
40  5.25  0.009
11  5.28  0.009
49  5.42  0.008
28  5.56  0.008
46  5.73  0.007
19  5.86  0.007
14  5.90  0.006
37  5.93  0.006
43  5.97  0.006
22  6.09  0.006
7   6.26  0.005
25  6.51  0.005
60  6.55  0.005
16  6.61  0.004
48  7.13  0.003
24  7.14  0.003
51  7.17  0.003
54  7.29  0.003
57  7.37  0.003
13  7.41  0.003
39  7.42  0.003
21  7.58  0.003
45  7.69  0.003
27  8.28  0.002
30  8.40  0.002
56  8.45  0.002
53  8.86  0.001
15  8.88  0.001
62  8.89  0.001
59  9.36  0.001
23  9.38  0.001
32  9.60  0.001
47  9.62  0.001
29 10.14  0.001
64 10.53  0.001
55 10.72  0.001
61 11.07  0.000
31 11.99  0.000
63 13.10  0.000
Models ranked by AICc(x) 
Code
# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2)) 

mr_averaged_aic2

Call:
model.avg(object = candidates, subset = delta < 2)

Component model call: 
rma.mv(yi = yi, V = VCV, mods = ~<7 unique rhs>, random = list(~1 | 
     effectID, ~1 | paperID, ~1 | species_cleaned), data = dat, method = ML, 
     test = t, sparse = TRUE)

Component models: 
       df  logLik   AICc delta weight
1       6 -276.83 565.78  0.00   0.24
12      7 -275.93 566.03  0.25   0.21
14      7 -276.28 566.72  0.94   0.15
124     8 -275.48 567.17  1.39   0.12
(Null)  4 -279.77 567.59  1.81   0.10
134     9 -274.70 567.67  1.89   0.09
13      8 -275.78 567.78  2.00   0.09

Term codes: 
     age_class_clean      dispersal_phase fitness_higher_level 
                   1                    2                    3 
       whose_fitness 
                   4 

Model-averaged coefficients:  
(full average) 
                                 Estimate Std. Error z value Pr(>|z|)
intrcpt                          -0.08622    0.11029   0.782    0.434
age_class_cleanjuvenile           0.06241    0.05321   1.173    0.241
age_class_cleanmix               -0.10467    0.07067   1.481    0.139
dispersal_phasesettlement         0.02130    0.04137   0.515    0.607
whose_fitnessoffspring            0.01788    0.03542   0.505    0.614
fitness_higher_levelreproduction  0.03936    0.11059   0.356    0.722
fitness_higher_levelsurvival      0.03285    0.10117   0.325    0.745
 
(conditional average) 
                                 Estimate Std. Error z value Pr(>|z|)  
intrcpt                          -0.08622    0.11029   0.782   0.4344  
age_class_cleanjuvenile           0.06911    0.05169   1.337   0.1812  
age_class_cleanmix               -0.11591    0.06502   1.783   0.0746 .
dispersal_phasesettlement         0.06429    0.04902   1.312   0.1896  
whose_fitnessoffspring            0.04929    0.04369   1.128   0.2592  
fitness_higher_levelreproduction  0.21658    0.17002   1.274   0.2027  
fitness_higher_levelsurvival      0.18075    0.17198   1.051   0.2933  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# relative importance of each predictor for all the models
importance <- sw(candidates)

importance
                     age_class_clean dispersal_phase whose_fitness
Sum of weights:      0.74            0.40            0.40         
N containing models:   32              32              32         
                     fitness_higher_level sex  dispersal_type
Sum of weights:      0.28                 0.20 0.14          
N containing models:   32                   32   32          

7 Publication bias

Code
# if each group n is available - assume n/2
dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n,  
                      (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))
  
dat$sqeffectN <- sqrt(dat$effectN)

mod_egger <- rma.mv(yi = yi, 
                    V = VCV,
                mod = ~ sqeffectN,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_egger)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.1824   486.3649   496.3649   519.2263   496.4495   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0734  0.2710    717     no         effectID 
sigma^2.2  0.0221  0.1487    210     no          paperID 
sigma^2.3  0.0155  0.1244    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 15756.7080, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 715) = 0.0164, p-val = 0.8982

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0272  0.0400  -0.6800  715  0.4968  -0.1057  0.0513    
sqeffectN    0.0005  0.0040   0.1280  715  0.8982  -0.0074  0.0084    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_egger)*100, 2)
   R2_marginal R2_conditional 
          0.00          33.88 
Code
small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)
Code
small

Figure S26. Relationship between sample size and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# creating publication year from "reference"
dat$reference <- as.character(dat$reference)
dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))
dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))
# decline effect
mod_dec <- rma.mv(yi = yi, V = vi,
                     mod = ~ year,
                     data = dat, 
                     random = list(
                       ~ 1 | effectID,
                       ~ 1 | paperID,
                       ~ 1 | species_cleaned),
                     test = "t",
                     sparse = TRUE)
summary(mod_dec)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-223.4816   446.9632   456.9632   479.8246   457.0479   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0605  0.2460    717     no         effectID 
sigma^2.2  0.0161  0.1267    210     no          paperID 
sigma^2.3  0.0204  0.1428    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 715) = 7980.4347, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 715) = 0.3932, p-val = 0.5308

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -1.8896  2.9797  -0.6341  715  0.5262  -7.7396  3.9605    
year       0.0009  0.0015   0.6271  715  0.5308  -0.0020  0.0038    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_dec)*100, 2)
   R2_marginal R2_conditional 
          0.14          37.68 
Code
decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)
Code
decline 

Figure S27. Relationship between publication year and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# All together

mod_comb <- rma.mv(yi = yi, 
                   V = VCV,
                   mod = ~ year + sqeffectN,
                   data = dat, 
                   random = list(
                     ~ 1 | effectID,
                     ~ 1 | paperID,
                     ~ 1 | species_cleaned),
                   # ~ 1 | phylogeny),
                   # R= list(phylogeny = cor_tree),
                   test = "t",
                   sparse = TRUE)
                   
summary(mod_comb)

Multivariate Meta-Analysis Model (k = 717; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.7662   485.5323   497.5323   524.9576   497.6512   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0731  0.2705    717     no         effectID 
sigma^2.2  0.0202  0.1422    210     no          paperID 
sigma^2.3  0.0190  0.1379    150     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 714) = 15746.7656, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 714) = 0.1606, p-val = 0.8517

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -1.7908  3.1712  -0.5647  714  0.5724  -8.0167  4.4351    
year         0.0009  0.0016   0.5559  714  0.5785  -0.0022  0.0040    
sqeffectN    0.0001  0.0041   0.0295  714  0.9765  -0.0079  0.0081    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_comb)*100, 2)
   R2_marginal R2_conditional 
          0.11          35.00 
Code
# combined sample size and year
bubble_plot(mod_comb,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S28. Combined effect of sample size and publication year on effect size (plotting sample size only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# combined sample size and year
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S29. Combined effect of sample size and publication year on effect size (plotting publication year only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
dat <- dat %>%
  mutate(leave_out = reference)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(model) {
  df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
  return(df)
}

# using dplyr to form data frame
MA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout$left_out <- as.factor(MA_oneout$left_out)
MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)

# saving the runs
saveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))
Code
MA_oneout <- readRDS(here("Rdata", "MA_oneout.RDS"))

# plotting
leaveoneout <- ggplot(MA_oneout) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Study left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
Code
leaveoneout

Figure S30. Leave one study out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals.
Code
dat <- dat %>%
  mutate(leave_out = species_cleaned)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# # writing function for extracting est, ci.lb, and ci.ub from all models
# est.func <- function(model) {
#   df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
#   return(df)
# }

# using dplyr to form data frame
MA_oneout1 <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout1$left_out <- as.factor(MA_oneout1$left_out)
MA_oneout1$left_out <- factor(MA_oneout1$left_out, levels = MA_oneout1$left_out)

# saving the runs
saveRDS(MA_oneout1, here("Rdata", "MA_oneout1.RDS"))
Code
MA_oneout1 <- readRDS(here("Rdata", "MA_oneout1.RDS"))

# plotting
leaveoneout1 <- ggplot(MA_oneout1) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Species left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
Code
leaveoneout1

Figure S31. Leave one species out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals.
Code
precision_plot<-funnel(mod1, 
       yaxis="seinv",
       # = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 16),
       xlim = c(-10,15),
       digits=c(0,1))

Figure S32. Relationship between standardized residuals and precision (inverse of SE).

8 R Session Information

Code
# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.2.3 (2023-03-15)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: cowplot(v.1.1.1), matrixcalc(v.1.0-6), orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.1.3), ggstance(v.0.3.6), ggtree(v.3.6.2), ape(v.5.7-1), metafor(v.4.4-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.5-3), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.8.9), MuMIn(v.1.47.5), kableExtra(v.1.3.4), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.4) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): nlme(v.3.1-162), fs(v.1.6.3), webshot(v.0.5.5), RColorBrewer(v.1.1-3), httr(v.1.4.7), rprojroot(v.2.0.4), latex2exp(v.0.9.6), tools(v.4.2.3), utf8(v.1.2.3), R6(v.2.5.1), vipor(v.0.4.7), mgcv(v.1.8-42), lazyeval(v.0.2.2), withr(v.2.5.2), tidyselect(v.1.2.0), compiler(v.4.2.3), cli(v.3.6.1), rvest(v.1.0.3), pacman(v.0.5.1), xml2(v.1.3.5), sandwich(v.3.1-0), labeling(v.0.4.3), scales(v.1.4.0), mvtnorm(v.1.2-3), systemfonts(v.1.0.5), digest(v.0.6.33), yulab.utils(v.0.1.1), rmarkdown(v.2.25), svglite(v.2.1.2), pkgconfig(v.2.0.3), htmltools(v.0.5.6.1), highr(v.0.10), fastmap(v.1.1.1), htmlwidgets(v.1.6.1), rlang(v.1.1.1), rstudioapi(v.0.15.0), gridGraphics(v.0.5-1), farver(v.2.1.1), generics(v.0.1.3), zoo(v.1.8-12), jsonlite(v.1.8.7), magrittr(v.2.0.3), ggplotify(v.0.1.2), ggbeeswarm(v.0.7.2), Rcpp(v.1.0.11), fansi(v.1.0.5), lifecycle(v.1.0.4), stringi(v.1.7.12), multcomp(v.1.4-25), yaml(v.2.3.7), mathjaxr(v.1.6-0), MASS(v.7.3-58.2), grid(v.4.2.3), parallel(v.4.2.3), lattice(v.0.20-45), splines(v.4.2.3), hms(v.1.1.3), knitr(v.1.45), pillar(v.1.9.0), estimability(v.1.4.1), codetools(v.0.2-19), stats4(v.4.2.3), glue(v.1.6.2), evaluate(v.0.23), ggfun(v.0.1.3), vctrs(v.0.6.4), treeio(v.1.22.0), tzdb(v.0.4.0), gtable(v.0.3.6), cachem(v.1.0.8), xfun(v.0.40), xtable(v.1.8-4), tidytree(v.0.4.5), coda(v.0.19-4), survival(v.3.5-3), viridisLite(v.0.4.2), aplot(v.0.2.2), beeswarm(v.0.4.0), memoise(v.2.0.1), timechange(v.0.2.0) and TH.data(v.1.1-2)